Incremento del nivel medio del mar

Juan Sebastián Mendoza Páez

Estudiante Ing. Forestal

¿Por qué es importante?

NASA

“Desde 1880, el nivel del mar global ha aumentado 20 centimetros. Para el 2100 se proyecta que aumente entre 30 y 122 centimetros más.”

Bases de datos

Fuentes

  • Oficina Nacional de Administración Oceánica y Atmosférica (NOAA)
  • Organización de Investigación Científica e Industrial del Commonwealth (CSIRO)
Mostrar Código
pacman::p_load(TSstudio, tsbox, tsoutliers, tidyverse,
               lubridate, xts, magrittr, forecast, dygraphs)

url <- 'https://gml.noaa.gov/webdata/ccgg/trends/co2/co2_mm_mlo.txt'
cov <- read.table(file = url) %>%
  select(3, 4) %>%
  rename(fecha = 1, co2 = 2) %>%
  mutate(fecha = date_decimal(fecha) %>%
           floor_date(unit = 'month')) %>%
  filter(fecha <= '2013-12-01')

data <- read.csv('gmsl.csv') %>%
  mutate(fecha = date_decimal(Time) %>%
           floor_date(unit = 'month')) %>%
  rename(gsml = 2) %>%
  select(4, 2)%>%
  filter(fecha > '1958-02-01') %>%
  left_join(y = cov) %>%
  {xts(x = .[,2:3], order.by = .$fecha)}

rm(cov)

dygraph(data, main = 'Nivel medio del mar global vs
        Concentración de CO2 en la atmósfera') %>%
  dySeries("co2", axis = 'y2', label = 'CO2') %>%
  dySeries('gsml', label = 'GSML') %>%
  dyAxis('y', label = 'Nivel del mar (mm)') %>%
  dyAxis('y2', label = 'CO2 (ppm)', independentTicks = T) %>%
  dyRangeSelector()
Mostrar Código
ts_split(ts.obj = ts_ts(data), sample.out = 12) %T>%
  {assign(x = "train", value = .$train, envir = .GlobalEnv)} %T>%
  {assign(x = "test", value = .$test, envir = .GlobalEnv)}

Estacionalidad y prueba de raíces unitarias

Mostrar Código
acf(train[,'gsml'], main = '')

Mostrar Código
tseries::adf.test(x = train[,'gsml'])

    Augmented Dickey-Fuller Test

data:  train[, "gsml"]
Dickey-Fuller = -2.7778, Lag order = 8, p-value = 0.2491
alternative hypothesis: stationary

Modelo

Mostrar Código
matriz_diseño <- model.matrix(~ -1 + train[,'co2'])
colnames(matriz_diseño) <- 'co2'

modelo <- Arima(y = train[,'gsml'], 
                order = c(2,1,3),
                xreg = matriz_diseño,
                include.drift = T)

modelo %>% lmtest::coeftest()

z test of coefficients:

       Estimate Std. Error  z value  Pr(>|z|)    
ar1    0.797070   0.054991  14.4947 < 2.2e-16 ***
ar2   -0.095533   0.054919  -1.7395   0.08194 .  
ma1   -0.339090   0.040237  -8.4274 < 2.2e-16 ***
ma2    0.208116   0.038717   5.3753 7.646e-08 ***
ma3   -0.693544   0.028246 -24.5535 < 2.2e-16 ***
drift  0.208518   0.041471   5.0281 4.955e-07 ***
co2   -0.171228   0.087716  -1.9521   0.05093 .  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Modelo

Mostrar Código
equatiomatic::extract_eq(modelo, wrap = T, use_coefs = T)

\[ \begin{alignat}{2} &\text{let}\quad &&y_{t} = \operatorname{y}_{\operatorname{0}} +0.21\operatorname{t} -0.17\operatorname{co2}_{\operatorname{t}} +\eta_{t} \\ &\text{where}\quad &&(1 -0.8\operatorname{B} +0.1\operatorname{B}^{\operatorname{2}} )\ (1 - \operatorname{B}) \eta_{t} \\ & &&= (1 -0.34\operatorname{B} +0.21\operatorname{B}^{\operatorname{2}} -0.69\operatorname{B}^{\operatorname{3}} )\ \varepsilon_{t} \\ &\text{where}\quad &&\varepsilon_{t} \sim{WN(0, \sigma^{2})} \end{alignat} \]

Supuestos

Mostrar Código
modelo %>%
  checkresiduals(theme = theme_bw(), test = F)

Mostrar Código
modelo %>%
  residuals() %>%
  shapiro.test()

    Shapiro-Wilk normality test

data:  .
W = 0.99861, p-value = 0.8921
Mostrar Código
modelo %>%
  residuals() %>%
  Box.test(type = 'Ljung-Box')

    Box-Ljung test

data:  .
X-squared = 0.029686, df = 1, p-value = 0.8632

Outliers

Mostrar Código
tso(y = train[,'gsml'], xreg = train[,'co2'])
Series:  
Regression with ARIMA(1,1,3) errors 

Coefficients:
         ar1      ma1     ma2      ma3     xreg
      0.7152  -0.2650  0.1905  -0.6741  -0.1222
s.e.  0.0439   0.0387  0.0333   0.0275   0.0851

sigma^2 = 3.074:  log likelihood = -1299.53
AIC=2611.06   AICc=2611.19   BIC=2637.99

No outliers were detected.

Predicciones

Mostrar Código
prediccion <- forecast(object = modelo, h = 12, level = 95, xreg = test[,'co2'])

serie <- prediccion %>%
  as_tibble() %>%
  ts(start = c(2013, 1), end = c(2013, 12), frequency = 12) %>%
  {cbind(data[, 'gsml'], .)}

colnames(serie) <- c('gsml', 'predicho', 'lwr', 'upr')

serie %>%
  dygraph() %>%
  dySeries('gsml', label = 'Observada') %>% 
  dySeries(c('lwr', 'predicho', 'upr'), label = 'Predicho') %>%
  dyRangeSelector()
Mostrar Código
accuracy(prediccion, test[,'gsml'])[,c(2,5)]
                 RMSE     MAPE
Training set 1.719463      Inf
Test set     9.180928 12.04332

Gracias